********************************************************************************
* Estimating Social Preferences and Kantian Morality in Strategic Interactions *
* Boris van Leeuwen & Ingela Alger *********************************************
* This version: September 2024 *************************************************
********************************************************************************

clear

* set the directory
cd ".../SocialPrefs_KantianMorality_ReplicationPackage"

********************************************************************************
***** Table 1: Average actions and beliefs for each game protocol **************
********************************************************************************
use "hm_all.dta", clear

replace x3=. if game_type!="SPD"
replace y3=. if game_type!="SPD"

eststo clear
quietly estpost tabstat t r p s x1 x2 x3 y1 y2 y3, by(game_no)
esttab using "Tab1_sumstat_games.tex",label cells("t r p s x1(fmt(2)) x2(fmt(2)) x3(fmt(2)) y1(fmt(2)) y2(fmt(2)) y3(fmt(2))") nomtitle nonumber replace noobs nonote varwidth(5) drop("Total" )

*totals by game protocol type
quietly estpost tabstat x1 x2 x3 y1 y2 y3, by(game_type)
esttab,label cells("x1(fmt(2)) x2(fmt(2)) x3(fmt(2)) y1(fmt(2)) y2(fmt(2)) y3(fmt(2))") nomtitle nonumber replace noobs nonote varwidth(5) drop("Total" )

********************************************************************************
*** Figure 4: Distributions of individual parameter estimates ******************
********************************************************************************
use "hm_estimates.dta", clear

local estm = "abc_rn"

twoway__histogram_gen a_`estm' if core==1, frequency start(-2) width(.1) gen(a3_h a3_x)
twoway__histogram_gen b_`estm' if core==1, frequency start(-2) width(.1) gen(b3_h b3_x)
twoway__histogram_gen c_`estm' if core==1, frequency start(-2) width(.1) gen(c3_h c3_x)

** for the fitted Gumbel distributions, the area on the curve is normalized by multiplying with 11.2 (112*0.1=11.2) to make sure that the area under the curve is equal to the area of the bars. The values for "a" and "b" below come from the fitted Gumbel distributions (see online appendix).

	local a=0.0797 
	local b=0.1422
twoway (bar a3_h a3_x, barwidth(.1)) ///
	(function 11.2*((1/`b')*exp(-((x-`a')/`b')-exp(-((x-`a')/`b')))), range(-2 2) lcolor(black) lpattern(dash) lwidth(medthin)) ///
	, xlabel(-2(.4)2) xtitle(estimated {&alpha}{subscript:i}) ///
	ytitle(Frequency) ylabel(0(10)35) ///
	scheme(s1mono) xsize(6) ysize(2) scale(1.8) legend(off)	
graph save Graph "hist_a_`estm'_core2.gph", replace

	local a=-0.3193 
	local b=0.3314
twoway (bar b3_h b3_x, barwidth(.1)) ///
	(function 11.2*((1/`b')*exp(-(((-x)-`a')/`b')-exp(-(((-x)-`a')/`b')))), range(-2 2) lcolor(black) lpattern(dash) lwidth(medthin)) ///
	, xlabel(-2(.4)2) xtitle(estimated {&beta}{subscript:i}) ///
	ytitle(Frequency) ylabel(0(10)35) ///
	scheme(s1mono) xsize(6) ysize(2) scale(1.8) legend(off)		
graph save Graph "hist_b_`estm'_core2.gph", replace

	local a=0.0629 
	local b=0.1126
twoway (bar c3_h c3_x, barwidth(.1)) ///
	(function 11.2*((1/`b')*exp(-((x-`a')/`b')-exp(-((x-`a')/`b')))), range(-2 2) lcolor(black) lpattern(dash) lwidth(medthin)) ///
	, xlabel(-2(.4)2) xtitle(estimated {&kappa}{subscript:i}) ///
	ytitle(Frequency) ylabel(0(10)35) ///
	scheme(s1mono) xsize(6) ysize(2) scale(1.8) legend(off)		
graph save Graph "hist_c_`estm'_core2.gph", replace

graph combine "hist_a_`estm'_core2.gph" ///
	"hist_b_`estm'_core2.gph" ///
	"hist_c_`estm'_core2.gph" ///
		, colfirst cols(1) scheme(s1mono) commonscheme xsize(6) ysize(6) scale(.5)
graph export "fig4_distr_core.pdf", as(pdf) replace

erase "hist_a_`estm'_core2.gph"
erase "hist_b_`estm'_core2.gph"
erase "hist_c_`estm'_core2.gph"

signrank a_`estm'=0 if core==1
signrank b_`estm'=0 if core==1
signrank c_`estm'=0 if core==1

********************************************************************************
*** Table 3: Individual parameter estimates ************************************
********************************************************************************
use "hm_estimates.dta", clear

eststo clear
quietly estpost summarize a_abc_rn b_abc_rn c_abc_rn if core==1, listwise detail
esttab using "Tab3_sum_abc_rn_trunc.tex", label cell((p50(fmt(%9.2f)) mean(fmt(%9.2f)) sd(fmt(%9.2f)) min(fmt(%9.2f)) max(fmt(%9.2f)))) nonumber nomtitle replace

********************************************************************************
*** Figure 5: Correlations between estimated preference parameters *************
********************************************************************************
use "hm_estimates.dta", clear

local estm = "abc_rn"

twoway (scatter b_`estm' a_`estm' if core==1, mcolor(gray%40) msymbol(circle) msize(small)) ///
	(lfit b_`estm' a_`estm' if core==1, lcolor(black) lpattern(vshortdash)) ///
	, ytitle(estimated {&beta}{subscript:i}) yline(0) xtitle(estimated {&alpha}{subscript:i}) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_a_b_`estm'.gph", replace

twoway (scatter c_`estm' a_`estm' if core==1, mcolor(gray%40) msymbol(circle) msize(small)) ///
	(lfit c_`estm' a_`estm' if core==1, lcolor(black) lpattern(vshortdash)) ///
	, ytitle(estimated {&kappa}{subscript:i}) yline(0) xtitle(estimated {&alpha}{subscript:i}) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_a_c_`estm'.gph", replace

twoway (scatter c_`estm' b_`estm' if core==1, mcolor(gray%40) msymbol(circle) msize(small)) ///
	(lfit c_`estm' b_`estm' if core==1, lcolor(black) lpattern(vshortdash)) ///
	, ytitle(estimated {&kappa}{subscript:i}) yline(0) xtitle(estimated {&beta}{subscript:i}) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_b_c_`estm'.gph", replace

graph combine "corr_a_b_`estm'.gph" ///
	"corr_a_c_`estm'.gph" ///
	"corr_b_c_`estm'.gph" ///
	, cols(3) scheme(s1mono) commonscheme xsize(6) ysize(2) scale(1.5)
graph export "fig5_corr.pdf", as(pdf) replace

erase "corr_a_b_`estm'.gph"
erase "corr_a_c_`estm'.gph"
erase "corr_b_c_`estm'.gph"

spearman a_`estm' b_`estm' if core==1
spearman a_`estm' c_`estm' if core==1
spearman b_`estm' c_`estm' if core==1

reg b_`estm' a_`estm' if core==1
reg c_`estm' a_`estm' if core==1
reg c_`estm' b_`estm' if core==1

* reported in Appendix A3 (Copula estimations)
gen negb_`estm'= -b_`estm'
ktau a_`estm' negb_`estm' if core==1
ktau a_`estm' c_`estm' if core==1
ktau negb_`estm' c_`estm' if core==1
drop negb_`estm'

********************************************************************************
*** Figure 6: Posterior probabilities of type classifications ******************
********************************************************************************
use "taus.dta", clear

* 2 types
* NOTE THAT TYPE 1 AND 2 ARE REVERSED IN THE PAPER AND THE DATA
histogram tau_2_2t_abc_rn, bin(10) start(0) fraction ///
	xlabel(0(.2)1) ///
	xtitle({&tau}{subscript:i}{subscript:1}) ylabel(0(.2)1) ///
	scheme(s1mono) xsize(4) ysize(4) scale(1)	aspect(1) title("Type 1")
graph save Graph "hist_tau_rn_2_1.gph", replace

histogram tau_1_2t_abc_rn, bin(10) start(0) fraction ///
	xlabel(0(.2)1) ///
	xtitle({&tau}{subscript:i}{subscript:2}) ylabel(0(.2)1) ///
	scheme(s1mono) xsize(4) ysize(4) scale(1)	aspect(1) title("Type 2")
graph save Graph "hist_tau_rn_2_2.gph", replace 

graph combine "hist_tau_rn_2_1.gph" ///
	"hist_tau_rn_2_2.gph" ///
	, colfirst rows(1) scheme(s1mono) commonscheme xsize(12) ysize(4) scale(1) title("2-types model")
graph save Graph "hist_tau_rn_2types.gph", replace

erase "hist_tau_rn_2_1.gph"
erase "hist_tau_rn_2_2.gph"

* 3 types
* NOTE THAT TYPE 1 AND 3 ARE REVERSED IN THE PAPER AND THE DATA
histogram tau_3_3t_abc_rn, bin(10) start(0) fraction ///
	xlabel(0(.2)1) ///
	xtitle({&tau}{subscript:i}{subscript:1}) ylabel(0(.2)1) ///
	scheme(s1mono) xsize(4) ysize(4) scale(1)	aspect(1) title("Type 1")
graph save Graph "hist_tau_rn_3_1.gph", replace

histogram tau_2_3t_abc_rn, bin(10) start(0) fraction ///
	xlabel(0(.2)1) ///
	xtitle({&tau}{subscript:i}{subscript:2}) ylabel(0(.2)1) ///
	scheme(s1mono) xsize(4) ysize(4) scale(1)	aspect(1) title("Type 2")
graph save Graph "hist_tau_rn_3_2.gph", replace

histogram tau_1_3t_abc_rn, bin(10) start(0) fraction ///
	xlabel(0(.2)1) ///
	xtitle({&tau}{subscript:i}{subscript:3}) ylabel(0(.2)1) ///
	scheme(s1mono) xsize(4) ysize(4) scale(1)	aspect(1) title("Type 3")
graph save Graph "hist_tau_rn_3_3.gph", replace

graph combine "hist_tau_rn_3_1.gph" ///
	"hist_tau_rn_3_2.gph" ///
	"hist_tau_rn_3_3.gph" ///
	, colfirst rows(1) scheme(s1mono) commonscheme xsize(12) ysize(4) title("3-types model")
graph save Graph "hist_tau_rn_3types.gph", replace

erase "hist_tau_rn_3_1.gph"
erase "hist_tau_rn_3_2.gph"
erase "hist_tau_rn_3_3.gph"

graph combine "hist_tau_rn_2types.gph" ///
	"hist_tau_rn_3types.gph" ///
	, colfirst rows(2) scheme(s1mono) commonscheme xsize(6) ysize(5) scale(1.2)
graph export "fig6_hist_tau.pdf", as(pdf) replace	

erase "hist_tau_rn_2types.gph"
erase "hist_tau_rn_3types.gph"

********************************************************************************
*** Figure 7: ICL scores *******************************************************
********************************************************************************
loc i = 1 
loc linecolors = `"lc("251 109 76" "194 59 34" "138 0 0" "72 0 0")"'
loc linecolors2 = `"lc("0 153 229" "23 63 95" "52 191 73" "255 76 76")"'

* The data below are the ICL scores of the respective models
loc ab1 = 4968.63912609344
loc ab2 = 4611.65200582451
loc ab3 = 4567.02013782731

loc abde1 = 4957.56702151085
loc abde2 = 4533.26403498156
loc abde3 = 4470.57098722338

loc abc1 = 4901.07624231124
loc abc2 = 4545.04051963268
loc abc3 = 4504.75133777813

loc abcde1 = 4895.25919133067
loc abcde2 = 4471.8613557445
loc abcde3 = 4417.16120178324

twoway (scatteri `ab1' 1 `ab2' 2 `ab3' 3, recast(connected)) (scatteri `abde1' 1 `abde2' 2 `abde3' 3, recast(connected)) (scatteri `abc1' 1 `abc2' 2 `abc3' 3, recast(connected)) (scatteri `abcde1' 1 `abcde2' 2 `abcde3' 3, recast(connected)),  scheme(s1mono)	xla(1 (1) 3, ang(h))  xtitle("Number of types") ytitle("ICL") ylabel(4300(200)5000) legend(order(1 "{&alpha}, {&beta}" 2 "{&alpha}, {&beta}, {&delta}, {&gamma}" 3 "{&alpha}, {&beta}, {&kappa}" 4 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(2) ring(0)  rowgap(1)) xsize(5) ysize(4) 
graph export "fig7_icl.pdf", as(pdf) replace

********************************************************************************
*** Figure 8: Distributions of individual parameter estimates (incl recipr) ****
********************************************************************************
use "hm_estimates.dta", clear

local estm = "abcde_rn"

gen a_squeezed = a_`estm'
replace a_squeezed = 2.3 if a_`estm' > 2
replace a_squeezed = -2.3 if a_`estm' < -2

gen b_squeezed = b_`estm'
replace b_squeezed = 2.3 if b_`estm' > 2
replace b_squeezed = -2.3 if b_`estm' < -2

gen c_squeezed = c_`estm'
replace c_squeezed = 2.3 if c_`estm' > 2
replace c_squeezed = -2.3 if c_`estm' < -2

gen d_squeezed = d_`estm'
replace d_squeezed = 2.3 if d_`estm' > 2
replace d_squeezed = -2.3 if d_`estm' < -2

gen e_squeezed = e_`estm'
replace e_squeezed = 2.3 if e_`estm' > 2
replace e_squeezed = -2.3 if e_`estm' < -2

twoway__histogram_gen a_squeezed if core==1, frequency start(-2.3) width(.1) gen(a_h a_x)
twoway__histogram_gen b_squeezed if core==1, frequency start(-2.3) width(.1) gen(b_h b_x)
twoway__histogram_gen c_squeezed if core==1, frequency start(-2.3) width(.1) gen(c_h c_x)
twoway__histogram_gen d_squeezed if core==1, frequency start(-2.3) width(.1) gen(d_h d_x)
twoway__histogram_gen e_squeezed if core==1, frequency start(-2.3) width(.1) gen(e_h e_x)

twoway (bar a_h a_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&alpha}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(6) ysize(2) scale(1.8)	
graph save Graph "hist_a_`estm'_sq.gph", replace

twoway (bar b_h b_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&beta}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(6) ysize(2) scale(1.8)	
graph save Graph "hist_b_`estm'_sq.gph", replace

twoway (bar c_h c_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&kappa}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(6) ysize(2) scale(1.8)	
graph save Graph "hist_c_`estm'_sq.gph", replace

twoway (bar d_h d_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&delta}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(6) ysize(2) scale(1.8)	
graph save Graph "hist_d_`estm'_sq.gph", replace

twoway (bar e_h e_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&gamma}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(6) ysize(2) scale(1.8)	
graph save Graph "hist_e_`estm'_sq.gph", replace

graph combine "hist_a_`estm'_sq.gph" ///
	"hist_b_`estm'_sq.gph" ///
	"hist_c_`estm'_sq.gph" ///
	"hist_d_`estm'_sq.gph" ///
	"hist_e_`estm'_sq.gph" ///
	, colfirst cols(1) scheme(s1mono) commonscheme xsize(6) ysize(8) scale(.5)
graph export "fig8_distr_abcde.pdf", as(pdf) replace

erase "hist_a_`estm'_sq.gph"
erase "hist_b_`estm'_sq.gph"
erase "hist_c_`estm'_sq.gph"
erase "hist_d_`estm'_sq.gph"
erase "hist_e_`estm'_sq.gph"

********************************************************************************
*** Table 7: Best individual fit ***********************************************
********************************************************************************
use "hm_estimates.dta", clear

gen min_aic_rn=min(aic_abcde_rn, aic_abcd_rn, aic_abce_rn, aic_abde_rn, aic_abc_rn, aic_abd_rn, aic_abe_rn, aic_acd_rn, aic_bce_rn, aic_ab_rn, aic_ac_rn, aic_ad_rn, aic_bc_rn, aic_be_rn, aic_a_rn, aic_b_rn, aic_c_rn, aic_si_rn)
gen best_aic_rn="";
replace best_aic_rn="{full model}" if aic_abcde_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\beta,\kappa,\gamma$}" if aic_abce_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\beta,\delta,\gamma$}" if aic_abde_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\beta,\gamma$}" if aic_abe_rn==min_aic_rn
replace best_aic_rn="{$\beta,\kappa,\gamma$}" if aic_bce_rn==min_aic_rn
replace best_aic_rn="{$\beta,\gamma$}" if aic_be_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\beta,\kappa,\delta$}" if aic_abcd_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\beta,\kappa$}" if aic_abc_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\beta,\delta$}" if aic_abd_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\kappa,\delta$}" if aic_acd_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\beta$}" if aic_ab_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\kappa$}" if aic_ac_rn==min_aic_rn
replace best_aic_rn="{$\alpha,\delta$}" if aic_ad_rn==min_aic_rn
replace best_aic_rn="{$\beta,\kappa$}" if aic_bc_rn==min_aic_rn
replace best_aic_rn="{$\alpha$}" if aic_a_rn==min_aic_rn
replace best_aic_rn="{$\beta$}" if aic_b_rn==min_aic_rn
replace best_aic_rn="{$\kappa$}" if aic_c_rn==min_aic_rn
replace best_aic_rn="{-}" if aic_si_rn==min_aic_rn

gen min_bic_rn=min(bic_abcde_rn, bic_abcd_rn, bic_abce_rn, bic_abde_rn, bic_abc_rn, bic_abd_rn, bic_abe_rn, bic_acd_rn, bic_bce_rn, bic_ab_rn, bic_ac_rn, bic_ad_rn, bic_bc_rn, bic_be_rn, bic_a_rn, bic_b_rn, bic_c_rn, bic_si_rn)
gen best_bic_rn="";
replace best_bic_rn="{full model}" if bic_abcde_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\beta,\kappa,\gamma$}" if bic_abce_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\beta,\delta,\gamma$}" if bic_abde_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\beta,\gamma$}" if bic_abe_rn==min_bic_rn
replace best_bic_rn="{$\beta,\kappa,\gamma$}" if bic_bce_rn==min_bic_rn
replace best_bic_rn="{$\beta,\gamma$}" if bic_be_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\beta,\kappa,\delta$}" if bic_abcd_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\beta,\kappa$}" if bic_abc_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\beta,\delta$}" if bic_abd_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\kappa,\delta$}" if bic_acd_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\beta$}" if bic_ab_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\kappa$}" if bic_ac_rn==min_bic_rn
replace best_bic_rn="{$\alpha,\delta$}" if bic_ad_rn==min_bic_rn
replace best_bic_rn="{$\beta,\kappa$}" if bic_bc_rn==min_bic_rn
replace best_bic_rn="{$\alpha$}" if bic_a_rn==min_bic_rn
replace best_bic_rn="{$\beta$}" if bic_b_rn==min_bic_rn
replace best_bic_rn="{$\kappa$}" if bic_c_rn==min_bic_rn
replace best_bic_rn="{-}" if bic_si_rn==min_bic_rn

eststo clear
estpost tabulate best_bic_rn if core==1
esttab using "Tab7_bic_model_comparison_rn.tex" ///
	, cells("b(label(freq)) pct(fmt(1))") varlabels(, blist(Total "{hline @width}{break}")) ///
	nonumber nomtitle noobs ///
	order( ///
"{full model}" ///
"{$\alpha,\beta,\kappa,\delta$}" ///
"{$\alpha,\beta,\kappa,\gamma$}" ///
"{$\alpha,\beta,\delta,\gamma$}" ///
"{$\alpha,\beta,\kappa$}" ///
"{$\alpha,\beta,\delta$}" ///
"{$\alpha,\beta,\gamma$}" ///
"{$\alpha,\kappa,\delta$}" ///
"{$\beta,\kappa,\gamma$}" ///
"{$\alpha,\beta$}" ///
"{$\alpha,\kappa$}" ///
"{$\alpha,\delta$}" ///
"{$\beta,\kappa$}" ///
"{$\beta,\gamma$}" ///
"{$\alpha$}" ///
"{$\beta$}" ///
"{$\kappa$}" ///
"{-}" ///
) varwidth(40) substitute(\_ _) replace

eststo clear
estpost tabulate best_aic_rn if core==1
esttab using "Tab7_aic_model_comparison_rn.tex" ///
	, cells("b(label(freq)) pct(fmt(1))") varlabels(, blist(Total "{hline @width}{break}")) ///
	nonumber nomtitle noobs ///
	order( ///
"{full model}" ///
"{$\alpha,\beta,\kappa,\delta$}" ///
"{$\alpha,\beta,\kappa,\gamma$}" ///
"{$\alpha,\beta,\delta,\gamma$}" ///
"{$\alpha,\beta,\kappa$}" ///
"{$\alpha,\beta,\delta$}" ///
"{$\alpha,\beta,\gamma$}" ///
"{$\alpha,\kappa,\delta$}" ///
"{$\beta,\kappa,\gamma$}" ///
"{$\alpha,\beta$}" ///
"{$\alpha,\kappa$}" ///
"{$\alpha,\delta$}" ///
"{$\beta,\kappa$}" ///
"{$\beta,\gamma$}" ///
"{$\alpha$}" ///
"{$\beta$}" ///
"{$\kappa$}" ///
"{-}" ///
) varwidth(40) substitute(\_ _) replace

********************************************************************************
*** Figure 9: Accuracy of out-of-sample predictions ****************************
********************************************************************************
use "predictions.dta", clear
merge m:1 id using "core.dta"

loc linepatterns = `"lpattern("solid" "solid" "shortdash" "shortdash" "dot")"'
loc linecolors = `"lc("gray" "0 0 0" "0 0 0" "gray" "0 0 0")"'
loc u = "game_no"

preserve
keep if core == 1
loc j1=2
loc j2=18
loc k1=.1
loc k2=.1
loc k3=.8

collapse choice_correct, by(model `u')

distplot choice_correct if (model=="si_rn" | model=="ab_rn" | model=="abde_rn" | model=="abc_rn" | model=="abcde_rn") ///
	, over(model) scheme(s1mono) c(J J J J J) ///
	yla(`j2' "`j2'" 0(`j1')`j2', ang(h)) ///
	xla(`k1'(`k2')`k3', ang(h)) `linepatterns' `linecolors' ///
	xtitle("correct choice") ///
	legend(order(5 "self-interest" 1 "{&alpha}, {&beta}" 4 "{&alpha}, {&beta}, {&delta}, {&gamma}" 2 "{&alpha}, {&beta}, {&kappa}" 3 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(10) ring(0)  rowgap(1)) ///
	frequency title("individual estimates")
graph save Graph "out_ch_ind_`u'.gph", replace

distplot choice_correct if (model=="si_3t_rn" | model=="ab_3t_rn" | model=="abde_3t_rn" | model=="abc_3t_rn" | model=="abcde_3t_rn") ///
	, over(model) scheme(s1mono) c(J J J J J) ///
	yla(`j2' "`j2'" 0(`j1')`j2', ang(h)) ///
	xla(`k1'(`k2')`k3', ang(h)) `linepatterns' `linecolors' ///
	xtitle("correct choice") ///
	legend(order(5 "self-interest" 1 "{&alpha}, {&beta}" 4 "{&alpha}, {&beta}, {&delta}, {&gamma}" 2 "{&alpha}, {&beta}, {&kappa}" 3 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(10) ring(0)  rowgap(1)) ///
	frequency title("3-types models")
graph save Graph "out_ch_3t_`u'.gph", replace

distplot choice_correct if (model=="si_2t_rn" | model=="ab_2t_rn" | model=="abde_2t_rn" | model=="abc_2t_rn" | model=="abcde_2t_rn") ///
	, over(model) scheme(s1mono) c(J J J J J) ///
	yla(`j2' "`j2'" 0(`j1')`j2', ang(h)) ///
	xla(`k1'(`k2')`k3', ang(h)) `linepatterns' `linecolors' ///
	xtitle("correct choice") ///
	legend(order(5 "self-interest" 1 "{&alpha}, {&beta}" 4 "{&alpha}, {&beta}, {&delta}, {&gamma}" 2 "{&alpha}, {&beta}, {&kappa}" 3 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(10) ring(0)  rowgap(1)) ///
	frequency title("2-types models")
graph save Graph "out_ch_2t_`u'.gph", replace

distplot choice_correct if (model=="si_1t_rn" | model=="ab_1t_rn" | model=="abde_1t_rn" | model=="abc_1t_rn" | model=="abcde_1t_rn") ///
	, over(model) scheme(s1mono) c(J J J J J) ///
	yla(`j2' "`j2'" 0(`j1')`j2', ang(h)) ///
	xla(`k1'(`k2')`k3', ang(h)) `linepatterns' `linecolors' ///
	xtitle("correct choice") ///
	legend(order(5 "self-interest" 1 "{&alpha}, {&beta}" 4 "{&alpha}, {&beta}, {&delta}, {&gamma}" 2 "{&alpha}, {&beta}, {&kappa}" 3 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(10) ring(0)  rowgap(1)) ///
	frequency title("1-type models")
graph save Graph "out_ch_1t_`u'.gph", replace

restore
 
graph combine	"out_ch_ind_`u'.gph"	///
				"out_ch_3t_`u'.gph" 	///
				"out_ch_2t_`u'.gph" 	///
				"out_ch_1t_`u'.gph" 	///
	, rows(2) scheme(s1mono) commonscheme xsize(6) ysize(4) scale(.9)
graph export "fig9_out.pdf", as(pdf) replace

erase "out_ch_ind_`u'.gph"
erase "out_ch_3t_`u'.gph"
erase "out_ch_2t_`u'.gph"
erase "out_ch_1t_`u'.gph"

********************************************************************************
***** Figure A.1: Correlation between actions en beliefs ***********************
********************************************************************************
use "hm_all.dta", clear

replace x3=. if game_type!="SPD"
replace y3=. if game_type!="SPD"

preserve
collapse x1 x2 x3 y1 y2 y3, by(game_no game_type)
reshape long x y, i(game_no) j(action)
twoway (scatter y x, mcolor("255 76 76 %40") msymbol(circle) msize(medium)) ///
	(scatteri 0 0 1 1, recast(line) lcolor(black) lpattern(vshortdash)) ///
	, ytitle("mean belief (y)") xtitle("mean action (x)") ///
	ylabel(0(.2)1) xlabel(0(.2)1) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph export "figA1_corr_actions_beliefs.pdf", as(pdf) replace
restore

********************************************************************************
***** Table A.1: Lottery choices ***********************************************
********************************************************************************
use "allsubjects.dta", clear
tab risk if session>1
tab risk if session==1

********************************************************************************
*** Figure A.2: Distributions of individual parameter estimates (all subjects) *
********************************************************************************
use "hm_estimates.dta", clear

local estm = "abc_rn"

gen a_squeezed = a_`estm'
replace a_squeezed = 2.3 if a_`estm' > 2
replace a_squeezed = -2.3 if a_`estm' < -2

gen b_squeezed = b_`estm'
replace b_squeezed = 2.3 if b_`estm' > 2
replace b_squeezed = -2.3 if b_`estm' < -2

gen c_squeezed = c_`estm'
replace c_squeezed = 2.3 if c_`estm' > 2
replace c_squeezed = -2.3 if c_`estm' < -2

twoway__histogram_gen a_squeezed, frequency start(-2.3) width(.1) gen(a_h a_x)
twoway__histogram_gen b_squeezed, frequency start(-2.3) width(.1) gen(b_h b_x)
twoway__histogram_gen c_squeezed, frequency start(-2.3) width(.1) gen(c_h c_x)

twoway (bar a_h a_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&alpha}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_a_`estm'_squeezed.gph", replace

twoway (bar b_h b_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&beta}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_b_`estm'_squeezed.gph", replace

twoway (bar c_h c_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&kappa}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_c_`estm'_squeezed.gph", replace

graph combine "hist_a_`estm'_squeezed.gph" ///
	"hist_b_`estm'_squeezed.gph" ///
	"hist_c_`estm'_squeezed.gph" ///
	, colfirst cols(1) scheme(s1mono) commonscheme xsize(6) ysize(6) scale(.5)
graph export "figA2_distr_all.pdf", as(pdf) replace

erase "hist_a_`estm'_squeezed.gph"
erase "hist_b_`estm'_squeezed.gph"
erase "hist_c_`estm'_squeezed.gph"

********************************************************************************
*** Table A.3: Individual parameter estimates (all subjects) *******************
********************************************************************************
use "hm_estimates.dta", clear

eststo clear
quietly estpost summarize a_abc_rn b_abc_rn c_abc_rn, listwise detail
esttab using "TabA3_sum_abc_rn.tex", label cell((p50(fmt(%9.2f)) mean(fmt(%9.2f)) sd(fmt(%9.2f)) min(fmt(%9.2f)) max(fmt(%9.2f)))) nonumber nomtitle replace

********************************************************************************
*** Table A.6: Transitions between types ***************************************
********************************************************************************
use "taus.dta", clear

*Panel A
* note that in the paper, the labeling of types 1 and 2 for the 2-types model are switched
* note that in the paper, the labeling of types 1 and 3 for the 3-types model are switched
tab type_3t_abc_rn type_2t_abc_rn

*Panel B
* note that in the paper, the labeling of types 1 and 2 for the 2-types risk-neutral model are switched
tab type_2t_abc_rn type_2t_abc_ln

*Panel C
* note that in the paper, the labeling of types 1 and 3 for the 3-types risk-neutral model are switched
tab type_3t_abc_rn type_3t_abc_ln

* Panel D
* note that in the paper, the labeling of types 1 and 2 for both 2-types models are switched 
tab type_2t_abc_rn type_2t_abc_rn_re

* Panel E
* note that in the paper, the labeling of types 1 and 3 for the 3-types risk-neutral subjective expectations model are switched
* note that in the paper, the labeling of the types in the RE model are switched (type 1 in data = type 2 in paper; type 2 in data = type 1 in paper; type 3 in data = type 3 in paper)
tab type_3t_abc_rn type_3t_abc_rn_re

* Panel F
* note that in the paper, the labeling of types 1 and 2 for the alpha,beta,kappa model are switched 
tab type_2t_abc_rn type_2t_abcde_rn

* Panel G
* note that in the paper, the labeling of types 1 and 3 for the 3-types alpha,beta,kappa model are switched
* note that in the paper, the labeling of types 1 and 3 for the 3-types alpha,beta,kappa,gamma,delta model are switched
tab type_3t_abc_rn type_3t_abcde_rn

* Panel H
tab type_2t_abde_rn type_2t_abcde_rn

* Panel I
* note that in the paper, the labeling of types 1 and 3 for the 3-types alpha,beta,kappa,gamma,delta model are switched
* note that in the paper, the labeling of types 1 and 3 for the 3-types alpha,beta,gamma,delta model are switched
tab type_3t_abde_rn type_3t_abcde_rn

********************************************************************************
*** Table A.7: Strategies by type (distributional pref., and Kantian morality) *
********************************************************************************
use "hm_all.dta", clear
joinby id using "taus.dta", unmatched(none)

local estm = "abc_rn"

label def spd_strategy 1 "C,C,C" 2 "C,C,D" 3 "C,D,C" 4 "C,D,D" 5 "D,C,C" 6 "D,C,D" 7 "D,D,C" 8 "D,D,D"
label def tg_strategy 1 "I,G" 2 "I,K" 3 "N,G" 4 "N,K" 
label def ug_strategy 1 "E,A" 2 "E,F" 3 "U,A" 4 "U,F" 

gen spd_choice = choice if game_type=="SPD"
gen tg_choice = choice if game_type=="TG"
gen ug_choice = choice if game_type=="UG"

label values spd_choice spd_strategy
label values tg_choice tg_strategy
label values ug_choice ug_strategy

* Note that for the 2-types model, Type 1 and Type 2 in the data and the paper are reversed
* Note that for the 3-types model, Type 1 and Type 3 in the data and the paper are reversed
eststo clear
quietly estpost tabulate spd_choice type_2t_`estm' if game_type == "SPD"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

eststo clear
quietly estpost tabulate spd_choice type_3t_`estm' if game_type == "SPD"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

eststo clear
quietly estpost tabulate tg_choice type_2t_`estm' if game_type == "TG"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

eststo clear
quietly estpost tabulate tg_choice type_3t_`estm' if game_type == "TG"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")
	
eststo clear
quietly estpost tabulate ug_choice type_2t_`estm' if game_type == "UG"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

eststo clear
quietly estpost tabulate ug_choice type_3t_`estm' if game_type == "UG"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

********************************************************************************
*** Table A.8: Strategies by type (distributional pref., Kantian morality, and reciprocity
********************************************************************************
use "hm_all.dta", clear
joinby id using "taus.dta", unmatched(none)

local estm = "abcde_rn"

label def spd_strategy 1 "C,C,C" 2 "C,C,D" 3 "C,D,C" 4 "C,D,D" 5 "D,C,C" 6 "D,C,D" 7 "D,D,C" 8 "D,D,D"
label def tg_strategy 1 "I,G" 2 "I,K" 3 "N,G" 4 "N,K" 
label def ug_strategy 1 "E,A" 2 "E,F" 3 "U,A" 4 "U,F" 

gen spd_choice = choice if game_type=="SPD"
gen tg_choice = choice if game_type=="TG"
gen ug_choice = choice if game_type=="UG"

label values spd_choice spd_strategy
label values tg_choice tg_strategy
label values ug_choice ug_strategy

* Note that in the paper, the labeling of types in the 3-types model are switched (type 1 in data = type 3 in paper; type 2 in data = type 2 in paper; type 3 in data = type 1 in paper)
eststo clear
quietly estpost tabulate spd_choice type_2t_`estm' if game_type == "SPD"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

eststo clear
quietly estpost tabulate spd_choice type_3t_`estm' if game_type == "SPD"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

eststo clear
quietly estpost tabulate tg_choice type_2t_`estm' if game_type == "TG"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

eststo clear
quietly estpost tabulate tg_choice type_3t_`estm' if game_type == "TG"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")
	
eststo clear
quietly estpost tabulate ug_choice type_2t_`estm' if game_type == "UG"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

eststo clear
quietly estpost tabulate ug_choice type_3t_`estm' if game_type == "UG"
esttab, cell(colpct(fmt(0))) unstack noobs nonumber nomtitle label  drop("Total")

********************************************************************************
*** Figure A.3: Distributions of individual parameter estimates (incl delta & gamma + all subjects) *
********************************************************************************
use "hm_estimates.dta", clear

local estm = "abcde_rn"

gen a_squeezed = a_`estm'
replace a_squeezed = 2.3 if a_`estm' > 2
replace a_squeezed = -2.3 if a_`estm' < -2

gen b_squeezed = b_`estm'
replace b_squeezed = 2.3 if b_`estm' > 2
replace b_squeezed = -2.3 if b_`estm' < -2

gen c_squeezed = c_`estm'
replace c_squeezed = 2.3 if c_`estm' > 2
replace c_squeezed = -2.3 if c_`estm' < -2

gen d_squeezed = d_`estm'
replace d_squeezed = 2.3 if d_`estm' > 2
replace d_squeezed = -2.3 if d_`estm' < -2

gen e_squeezed = e_`estm'
replace e_squeezed = 2.3 if e_`estm' > 2
replace e_squeezed = -2.3 if e_`estm' < -2

twoway__histogram_gen a_squeezed, frequency start(-2.3) width(.1) gen(a_h a_x)
twoway__histogram_gen b_squeezed, frequency start(-2.3) width(.1) gen(b_h b_x)
twoway__histogram_gen c_squeezed, frequency start(-2.3) width(.1) gen(c_h c_x)
twoway__histogram_gen d_squeezed, frequency start(-2.3) width(.1) gen(d_h d_x)
twoway__histogram_gen e_squeezed, frequency start(-2.3) width(.1) gen(e_h e_x)

twoway (bar a_h a_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&alpha}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_a_`estm'_squeezed.gph", replace

twoway (bar b_h b_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&beta}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_b_`estm'_squeezed.gph", replace

twoway (bar c_h c_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&kappa}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_c_`estm'_squeezed.gph", replace

twoway (bar e_h e_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&gamma}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_e_`estm'_squeezed.gph", replace

twoway (bar d_h d_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&delta}{subscript:i}) ylabel(0(10)40) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_d_`estm'_squeezed.gph", replace

graph combine "hist_a_`estm'_squeezed.gph" ///
	"hist_b_`estm'_squeezed.gph" ///
	"hist_c_`estm'_squeezed.gph" ///
	"hist_d_`estm'_squeezed.gph" ///
	"hist_e_`estm'_squeezed.gph" ///
	, colfirst cols(1) scheme(s1mono) commonscheme xsize(6) ysize(8) scale(.5)
graph export "figA3_distr_abcde_all.pdf", as(pdf) replace

erase "hist_a_`estm'_squeezed.gph"
erase "hist_b_`estm'_squeezed.gph"
erase "hist_c_`estm'_squeezed.gph"
erase "hist_d_`estm'_squeezed.gph"
erase "hist_e_`estm'_squeezed.gph"

********************************************************************************
*** Table A.12: Individual parameter estimates (incl recipr) *******************
********************************************************************************
use "hm_estimates.dta", clear

eststo clear
quietly estpost summarize a_abcde_rn b_abcde_rn c_abcde_rn d_abcde_rn e_abcde_rn if core==1, listwise detail
esttab using "TabA12_sum_abcde_rn_trunc.tex", label cell((p50(fmt(%9.2f)) mean(fmt(%9.2f)) sd(fmt(%9.2f)) min(fmt(%9.2f)) max(fmt(%9.2f)))) nonumber nomtitle replace

********************************************************************************
***** Table A.13: Average predictive accuracy of out-of-sample predictions *****
********************************************************************************
use "predictions.dta", clear
merge m:1 id using "core.dta"

table ( model ) () () if core==1 & strpos(model, "_rn"), statistic(mean choice_correct) nformat(%4.3f mean) showcounts

********************************************************************************
*** Figure A.4: Distributions of individual parameter estimates (risk aversion) *
********************************************************************************
use "hm_estimates.dta", clear

local estm = "abc_crra"

gen a_squeezed = a_`estm'
replace a_squeezed = 2.3 if a_`estm' > 2
replace a_squeezed = -2.3 if a_`estm' < -2

gen b_squeezed = b_`estm'
replace b_squeezed = 2.3 if b_`estm' > 2
replace b_squeezed = -2.3 if b_`estm' < -2

gen c_squeezed = c_`estm'
replace c_squeezed = 2.3 if c_`estm' > 2
replace c_squeezed = -2.3 if c_`estm' < -2

twoway__histogram_gen a_squeezed if core==1, frequency start(-2.3) width(.1) gen(a_h a_x)
twoway__histogram_gen b_squeezed if core==1, frequency start(-2.3) width(.1) gen(b_h b_x)
twoway__histogram_gen c_squeezed if core==1, frequency start(-2.3) width(.1) gen(c_h c_x)

twoway (bar a_h a_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&alpha}{subscript:i}) ylabel(0(10)30) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_a_`estm'_sq.gph", replace

twoway (bar b_h b_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&beta}{subscript:i}) ylabel(0(10)30) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_b_`estm'_sq.gph", replace

twoway (bar c_h c_x, barwidth(.1)), xlabel(-2(.4)2 -2.25 "<" 2.25 ">") ///
	xtitle(estimated {&kappa}{subscript:i}) ylabel(0(10)30) ///
	scheme(s1mono) xsize(12) ysize(4) scale(1.8)	
graph save Graph "hist_c_`estm'_sq.gph", replace

graph combine "hist_a_`estm'_sq.gph" ///
	"hist_b_`estm'_sq.gph" ///
	"hist_c_`estm'_sq.gph" ///
	, colfirst cols(1) scheme(s1mono) commonscheme xsize(6) ysize(6) scale(.5)
graph export "figA4_distr_abc_crra.pdf", as(pdf) replace

erase "hist_a_`estm'_sq.gph"
erase "hist_b_`estm'_sq.gph"
erase "hist_c_`estm'_sq.gph"

signrank a_abc_crra=0 if core==1
signrank b_abc_crra=0 if core==1
signrank c_abc_crra=0 if core==1

********************************************************************************
*** Figure A.5: Correlations between risk neutral and CRRA estimates ***********
********************************************************************************
use "hm_estimates.dta", clear

twoway (scatter a_abc_crra a_abc_rn, mcolor(midblue%40) msymbol(circle) msize(small)) ///
	(scatteri -2 -2 2 2, recast(line) lcolor(black) lpattern(vshortdash)) ///
	if core==1 & inrange(a_abc_crra, -2, 2) & inrange(a_abc_rn, -2, 2) ///
	, ytitle(estimated {&alpha}{subscript:i} (crra)) yline(0) xtitle(estimated {&alpha}{subscript:i} (risk-neutral)) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_a_rn_crra_squeezed.gph", replace

twoway (scatter b_abc_crra b_abc_rn, mcolor(midblue%40) msymbol(circle) msize(small)) ///
	(scatteri -2 -2 2 2, recast(line) lcolor(black) lpattern(vshortdash)) ///
	if core==1 & inrange(b_abc_crra, -2, 2) & inrange(b_abc_rn, -2, 2) ///
	, ytitle(estimated {&beta}{subscript:i} (crra)) yline(0) xtitle(estimated {&beta}{subscript:i} (risk-neutral)) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_b_rn_crra_squeezed.gph", replace

twoway (scatter c_abc_crra c_abc_rn, mcolor(midblue%40) msymbol(circle) msize(small)) ///
	(scatteri -2 -2 2 2, recast(line) lcolor(black) lpattern(vshortdash)) ///
	if core==1 & inrange(c_abc_crra, -2, 2) & inrange(c_abc_rn, -2, 2) ///
	, ytitle(estimated {&kappa}{subscript:i} (crra)) yline(0) xtitle(estimated {&kappa}{subscript:i} (risk-neutral)) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_c_rn_crra_squeezed.gph", replace

graph combine "corr_a_rn_crra_squeezed.gph" ///
	"corr_b_rn_crra_squeezed.gph" ///
	"corr_c_rn_crra_squeezed.gph" ///
	, colfirst cols(3) scheme(s1mono) commonscheme xsize(12) ysize(4) scale(1.8)
graph export "figA5_corr_rn_crra.pdf", as(pdf) replace

erase "corr_a_rn_crra_squeezed.gph"
erase "corr_b_rn_crra_squeezed.gph"
erase "corr_c_rn_crra_squeezed.gph"

spearman a_abc_crra a_abc_rn if core==1
spearman b_abc_crra b_abc_rn if core==1
spearman c_abc_crra c_abc_rn if core==1

signrank a_abc_crra=a_abc_rn if core==1
signrank b_abc_crra=b_abc_rn if core==1
signrank c_abc_crra=c_abc_rn if core==1

gen ll_larger_crra = 0
replace ll_larger_crra = 1 if ll_abc_crra>ll_abc_rn 
tab ll_larger_crra if core==1

********************************************************************************
*** Figure A.6: ICL scores (under CRRA) ****************************************
********************************************************************************

loc i = 1 
loc linecolors = `"lc("251 109 76" "194 59 34" "138 0 0" "72 0 0")"'
loc linecolors2 = `"lc("0 153 229" "23 63 95" "52 191 73" "255 76 76")"'

loc ab1 = 4907.27498497117
loc ab2 = 4538.96326944463
loc ab3 = 4492.94901228689

loc abde1 = 4884.42495313794
loc abde2 = 4471.64483568485
loc abde3 = 4420.30837828012

loc abc1 = 4732.45662629125
loc abc2 = 4368.21453950349
loc abc3 = 4329.76894084824

loc abcde1 = 4727.78573959934
loc abcde2 = 4340.62037106791
loc abcde3 = 4304.57644927769

twoway (scatteri `ab1' 1 `ab2' 2 `ab3' 3, recast(connected)) (scatteri `abde1' 1 `abde2' 2 `abde3' 3, recast(connected)) (scatteri `abc1' 1 `abc2' 2 `abc3' 3, recast(connected)) (scatteri `abcde1' 1 `abcde2' 2 `abcde3' 3, recast(connected)), ylabel(4300(200)5000) title("logarithmic utility") scheme(s1mono)	xla(1 (1) 3, ang(h))  xtitle("Number of types") ytitle("ICL") legend(order(1 "{&alpha}, {&beta}" 2 "{&alpha}, {&beta}, {&delta}, {&gamma}" 3 "{&alpha}, {&beta}, {&kappa}" 4 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(2) ring(0)  rowgap(1))
graph export "figA6_icl_ln.pdf", as(pdf) replace

********************************************************************************
*** Table A.16: Best individual fit (under CRRA) *******************************
********************************************************************************
use "hm_estimates.dta", clear

gen min_aic_crra=min(aic_abcde_crra, aic_abcd_crra, aic_abce_crra, aic_abde_crra, aic_abc_crra, aic_abd_crra, aic_abe_crra, aic_acd_crra, aic_bce_crra, aic_ab_crra, aic_ac_crra, aic_ad_crra, aic_bc_crra, aic_be_crra, aic_a_crra, aic_b_crra, aic_c_crra, aic_si_crra)
gen best_aic_crra="";
replace best_aic_crra="{full model}" if aic_abcde_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\beta,\kappa,\gamma$}" if aic_abce_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\beta,\delta,\gamma$}" if aic_abde_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\beta,\gamma$}" if aic_abe_crra==min_aic_crra
replace best_aic_crra="{$\beta,\kappa,\gamma$}" if aic_bce_crra==min_aic_crra
replace best_aic_crra="{$\beta,\gamma$}" if aic_be_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\beta,\kappa,\delta$}" if aic_abcd_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\beta,\kappa$}" if aic_abc_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\beta,\delta$}" if aic_abd_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\kappa,\delta$}" if aic_acd_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\beta$}" if aic_ab_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\kappa$}" if aic_ac_crra==min_aic_crra
replace best_aic_crra="{$\alpha,\delta$}" if aic_ad_crra==min_aic_crra
replace best_aic_crra="{$\beta,\kappa$}" if aic_bc_crra==min_aic_crra
replace best_aic_crra="{$\alpha$}" if aic_a_crra==min_aic_crra
replace best_aic_crra="{$\beta$}" if aic_b_crra==min_aic_crra
replace best_aic_crra="{$\kappa$}" if aic_c_crra==min_aic_crra
replace best_aic_crra="{-}" if aic_si_crra==min_aic_crra


gen min_bic_crra=min(bic_abcde_crra, bic_abcd_crra, bic_abce_crra, bic_abde_crra, bic_abc_crra, bic_abd_crra, bic_abe_crra, bic_acd_crra, bic_bce_crra, bic_ab_crra, bic_ac_crra, bic_ad_crra, bic_bc_crra, bic_be_crra, bic_a_crra, bic_b_crra, bic_c_crra, bic_si_crra)
gen best_bic_crra="";
replace best_bic_crra="{full model}" if bic_abcde_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\beta,\kappa,\gamma$}" if bic_abce_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\beta,\delta,\gamma$}" if bic_abde_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\beta,\gamma$}" if bic_abe_crra==min_bic_crra
replace best_bic_crra="{$\beta,\kappa,\gamma$}" if bic_bce_crra==min_bic_crra
replace best_bic_crra="{$\beta,\gamma$}" if bic_be_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\beta,\kappa,\delta$}" if bic_abcd_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\beta,\kappa$}" if bic_abc_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\beta,\delta$}" if bic_abd_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\kappa,\delta$}" if bic_acd_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\beta$}" if bic_ab_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\kappa$}" if bic_ac_crra==min_bic_crra
replace best_bic_crra="{$\alpha,\delta$}" if bic_ad_crra==min_bic_crra
replace best_bic_crra="{$\beta,\kappa$}" if bic_bc_crra==min_bic_crra
replace best_bic_crra="{$\alpha$}" if bic_a_crra==min_bic_crra
replace best_bic_crra="{$\beta$}" if bic_b_crra==min_bic_crra
replace best_bic_crra="{$\kappa$}" if bic_c_crra==min_bic_crra
replace best_bic_crra="{-}" if bic_si_crra==min_bic_crra

eststo clear
estpost tabulate best_bic_crra if core==1
esttab using "TabA16_bic_model_comparison_crra.tex" ///
	, cells("b(label(freq)) pct(fmt(1))") varlabels(, blist(Total "{hline @width}{break}")) ///
	nonumber nomtitle noobs ///
	order( ///
"{full model}" ///
"{$\alpha,\beta,\kappa,\delta$}" ///
"{$\alpha,\beta,\kappa,\gamma$}" ///
"{$\alpha,\beta,\delta,\gamma$}" ///
"{$\alpha,\beta,\kappa$}" ///
"{$\alpha,\beta,\delta$}" ///
"{$\alpha,\beta,\gamma$}" ///
"{$\alpha,\kappa,\delta$}" ///
"{$\beta,\kappa,\gamma$}" ///
"{$\alpha,\beta$}" ///
"{$\alpha,\kappa$}" ///
"{$\alpha,\delta$}" ///
"{$\beta,\kappa$}" ///
"{$\beta,\gamma$}" ///
"{$\alpha$}" ///
"{$\beta$}" ///
"{$\kappa$}" ///
"{-}" ///
) varwidth(40) substitute(\_ _) replace

eststo clear
estpost tabulate best_aic_crra if core==1
esttab using "TabA16_aic_model_comparison_crra.tex" ///
	, cells("b(label(freq)) pct(fmt(1))") varlabels(, blist(Total "{hline @width}{break}")) ///
	nonumber nomtitle noobs ///
	order( ///
"{full model}" ///
"{$\alpha,\beta,\kappa,\delta$}" ///
"{$\alpha,\beta,\kappa,\gamma$}" ///
"{$\alpha,\beta,\delta,\gamma$}" ///
"{$\alpha,\beta,\kappa$}" ///
"{$\alpha,\beta,\delta$}" ///
"{$\alpha,\beta,\gamma$}" ///
"{$\alpha,\kappa,\delta$}" ///
"{$\beta,\kappa,\gamma$}" ///
"{$\alpha,\beta$}" ///
"{$\alpha,\kappa$}" ///
"{$\alpha,\delta$}" ///
"{$\beta,\kappa$}" ///
"{$\beta,\gamma$}" ///
"{$\alpha$}" ///
"{$\beta$}" ///
"{$\kappa$}" ///
"{-}" ///
) varwidth(40) substitute(\_ _) replace

********************************************************************************
*** Figure A.7: Accuracy of out-of-sample predictions (under CRRA) *************
********************************************************************************
use "predictions.dta", clear
merge m:1 id using "core.dta"
 
loc linepatterns = `"lpattern("solid" "solid" "shortdash" "shortdash" "dot")"'
loc linecolors = `"lc("gray" "0 0 0" "0 0 0" "gray" "0 0 0")"'
loc u = "game_no"

preserve
keep if core == 1
loc j1=2
loc j2=18
loc k1=.1
loc k2=.1
loc k3=.8

collapse choice_correct, by(model `u')

distplot choice_correct if (model=="si_crra" | model=="ab_crra" | model=="abde_crra" | model=="abc_crra" | model=="abcde_crra") ///
	, over(model) scheme(s1mono) c(J J J J J) ///
	yla(`j2' "`j2'" 0(`j1')`j2', ang(h)) ///
	xla(`k1'(`k2')`k3', ang(h)) `linepatterns' `linecolors' ///
	xtitle("correct choice") ///
	legend(order(5 "self-interest" 1 "{&alpha}, {&beta}" 4 "{&alpha}, {&beta}, {&delta}, {&gamma}" 2 "{&alpha}, {&beta}, {&kappa}" 3 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(10) ring(0)  rowgap(1)) ///
	frequency title("individual estimates")
graph save Graph "out_ch_ind_`u'_crra.gph", replace

distplot choice_correct if (model=="si_3t_ln" | model=="ab_3t_ln" | model=="abde_3t_ln" | model=="abc_3t_ln" | model=="abcde_3t_ln") ///
	, over(model) scheme(s1mono) c(J J J J J) ///
	yla(`j2' "`j2'" 0(`j1')`j2', ang(h)) ///
	xla(`k1'(`k2')`k3', ang(h)) `linepatterns' `linecolors' ///
	xtitle("correct choice") ///
	legend(order(5 "self-interest" 1 "{&alpha}, {&beta}" 4 "{&alpha}, {&beta}, {&delta}, {&gamma}" 2 "{&alpha}, {&beta}, {&kappa}" 3 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(10) ring(0)  rowgap(1)) ///
	frequency title("3-types models")
graph save Graph "out_ch_3t_`u'_crra.gph", replace

distplot choice_correct if (model=="si_2t_ln" | model=="ab_2t_ln" | model=="abde_2t_ln" | model=="abc_2t_ln" | model=="abcde_2t_ln") ///
	, over(model) scheme(s1mono) c(J J J J J) ///
	yla(`j2' "`j2'" 0(`j1')`j2', ang(h)) ///
	xla(`k1'(`k2')`k3', ang(h)) `linepatterns' `linecolors' ///
	xtitle("correct choice") ///
	legend(order(5 "self-interest" 1 "{&alpha}, {&beta}" 4 "{&alpha}, {&beta}, {&delta}, {&gamma}" 2 "{&alpha}, {&beta}, {&kappa}" 3 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(10) ring(0)  rowgap(1)) ///
	frequency title("2-types models")
graph save Graph "out_ch_2t_`u'_crra.gph", replace

distplot choice_correct if (model=="si_1t_ln" | model=="ab_1t_ln" | model=="abde_1t_ln" | model=="abc_1t_ln" | model=="abcde_1t_ln") ///
	, over(model) scheme(s1mono) c(J J J J J) ///
	yla(`j2' "`j2'" 0(`j1')`j2', ang(h)) ///
	xla(`k1'(`k2')`k3', ang(h)) `linepatterns' `linecolors' ///
	xtitle("correct choice") ///
	legend(order(5 "self-interest" 1 "{&alpha}, {&beta}" 4 "{&alpha}, {&beta}, {&delta}, {&gamma}" 2 "{&alpha}, {&beta}, {&kappa}" 3 "{&alpha}, {&beta}, {&kappa}, {&delta}, {&gamma}") col(1) pos(10) ring(0)  rowgap(1)) ///
	frequency title("1-type models")
graph save Graph "out_ch_1t_`u'_crra.gph", replace

restore

graph combine	"out_ch_ind_`u'_crra.gph"	///
				"out_ch_3t_`u'_crra.gph" 	///
				"out_ch_2t_`u'_crra.gph" 	///
				"out_ch_1t_`u'_crra.gph" 	///
	, rows(2) scheme(s1mono) commonscheme xsize(5) ysize(3) scale(.9)
graph export "figA7_out_crra.pdf", as(pdf) replace

erase "out_ch_ind_`u'_crra.gph"
erase "out_ch_3t_`u'_crra.gph"
erase "out_ch_2t_`u'_crra.gph"
erase "out_ch_1t_`u'_crra.gph"

********************************************************************************
*** Table A.17 Average predictive accuracy of out-of-sample predictions (under CRRA) *
********************************************************************************
use "predictions.dta", clear
merge m:1 id using "core.dta"

table ( model ) () () if core==1 & (strpos(model, "_ln") | strpos(model, "_crra")), statistic(mean choice_correct) nformat(%4.3f mean) showcounts

********************************************************************************
*** Figure A.8: Correlations between estimates using subjective and rational expectations *
********************************************************************************
use "hm_estimates.dta", clear

twoway (scatter a_abc_rn a_abc_rn_re, mcolor(midblue%40) msymbol(circle) msize(small)) ///
	(scatteri -2 -2 2 2, recast(line) lcolor(black) lpattern(vshortdash)) ///
	if core==1 & inrange(a_abc_rn, -2, 2) & inrange(a_abc_rn_re, -2, 2) ///
	, ytitle(estimated {&alpha}{subscript:i} (subjective)) yline(0) xtitle(estimated {&alpha}{subscript:i} (re)) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_a_re_subj_squeezed.gph", replace

twoway (scatter b_abc_rn b_abc_rn_re, mcolor(midblue%40) msymbol(circle) msize(small)) ///
	(scatteri -2 -2 2 2, recast(line) lcolor(black) lpattern(vshortdash)) ///
	if core==1 & inrange(b_abc_rn, -2, 2) & inrange(b_abc_rn_re, -2, 2) ///
	, ytitle(estimated {&beta}{subscript:i} (subjective)) yline(0) xtitle(estimated {&beta}{subscript:i} (re)) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_b_re_subj_squeezed.gph", replace

twoway (scatter c_abc_rn c_abc_rn_re, mcolor(midblue%40) msymbol(circle) msize(small)) ///
	(scatteri -2 -2 2 2, recast(line) lcolor(black) lpattern(vshortdash)) ///
	if core==1 & inrange(c_abc_rn, -2, 2) & inrange(c_abc_rn_re, -2, 2) ///
	, ytitle(estimated {&kappa}{subscript:i} (subjective)) yline(0) xtitle(estimated {&kappa}{subscript:i} (re)) ///
	ylabel(-2(.5)2) xlabel(-2(.5)2 ) xline(0) scheme(s1mono) ///
	xsize(4) ysize(4) scale(1.2) aspect(1) legend(off)
graph save Graph "corr_c_re_subj_squeezed.gph", replace

graph combine "corr_a_re_subj_squeezed.gph" ///
	"corr_b_re_subj_squeezed.gph" ///
	"corr_c_re_subj_squeezed.gph" ///
	, colfirst cols(3) scheme(s1mono) commonscheme xsize(6) ysize(2) scale(1.8)
graph export "figA8_corr_re_subj.pdf", as(pdf) replace

erase "corr_a_re_subj_squeezed.gph"
erase "corr_b_re_subj_squeezed.gph"
erase "corr_c_re_subj_squeezed.gph"

spearman a_abc_rn a_abc_rn_re if core==1
spearman b_abc_rn b_abc_rn_re if core==1
spearman c_abc_rn c_abc_rn_re if core==1
